home *** CD-ROM | disk | FTP | other *** search
/ Aminet 2 / Aminet AMIGA CDROM (1994)(Walnut Creek)[Feb 1994][W.O. 44790-1].iso / Aminet / util / misc / Fudgit233.lha / Source / src / spline.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-14  |  2.4 KB  |  119 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #ifndef NOUNISTD_H
  4. #include <unistd.h>
  5. #endif
  6. #include "fudgit.h"
  7. #include "dalloca.h"
  8. #include "head.h"
  9.  
  10. static double *Xs, *Ys, *Y2s;
  11. static int SplData = 0;
  12. static int initspline(double *x, double *y, int n);
  13.  
  14.  
  15. int Ft_spline(double *x, double *y, double yp1, double ypn, int n)
  16. {
  17.     int i, k;
  18.     double p, qn, sig, un, *u;
  19.     extern double *Y2s;
  20.  
  21.     /* initialize vectors Xs, Ys and Y2s */
  22.     if (initspline(x, y, n) == ERRR) {
  23.         fputs("spline: Allocation error.\n", stderr);
  24.         return(ERRR);
  25.     }
  26.     ADVECTOR(u, 1, n-1, "spline", return);
  27.     if (yp1 > 0.99e30) {
  28.         Y2s[1] = u[1] = 0.0;
  29.     }
  30.     else {
  31.         Y2s[1] = -0.5;
  32.         u[1] = (3.0/(x[2]-x[1])) * ((y[2]-y[1])/(x[2]-x[1]) - yp1);
  33.     }
  34.     for (i=2; i<=n-1; i++) {
  35.         sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
  36.         p = sig*Y2s[i-1]+2.0;
  37.         Y2s[i] = (sig-1.0)/p;
  38.         u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
  39.         u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1]) / p;
  40.     }
  41.     if (ypn > 0.99e30) {
  42.         qn = un = 0.0;
  43.     }
  44.     else {
  45.         qn = 0.5;
  46.         un = (3.0/(x[n] - x[n-1])) * (ypn - (y[n] - y[n-1])/(x[n]-x[n-1]));
  47.     }
  48.     Y2s[n] = (un - qn * u[n-1]) / (qn * Y2s[n-1] + 1.0);
  49.     for (k=n-1; k>=1; k--) {
  50.         Y2s[k] = Y2s[k] * Y2s[k+1] + u[k];
  51.     }
  52.     SplData = n;
  53.     return(0);
  54. }
  55.  
  56. /* routine called from C-calculator mode */
  57.  
  58. double Ft_interp(double x)
  59. {
  60.     int klo, khi, k;
  61.     double h, b, a, y;
  62.     extern double *Xs, *Ys, *Y2s;
  63.     extern int SplData;
  64.  
  65.     if (!SplData) {
  66.         fputs("Math error: interp: Routine not initialized.\n", stderr);
  67.         Ft_catcher(ERRR);
  68.     }
  69.     klo = 1;
  70.     khi = SplData;
  71.     while (khi-klo > 1) {
  72.         k = (khi+klo) >> 1;
  73.         if (Xs[k] > x)
  74.             khi = k;
  75.         else
  76.             klo = k;
  77.     }
  78.     h = Xs[khi] - Xs[klo];
  79.     if (h == 0.0) {
  80.         fputs("Math error: interp: Routine not properly initialized.\n",
  81.         stderr);
  82.         Ft_catcher(ERRR);
  83.     }
  84.     a = (Xs[khi]-x)/h;
  85.     b = (x-Xs[klo])/h;
  86.     y = a * Ys[klo] + b * Ys[khi]
  87.     + (a*(a*a-1.0) * Y2s[klo] + b*(b*b-1.0) * Y2s[khi]) * (h*h) / 6.0;
  88.     return(y);
  89. }
  90.  
  91. static int initspline(double *x, double *y, int n)
  92. {
  93.     int i;
  94.     extern double *Xs, *Ys, *Y2s;
  95.     extern void Ft_free_dvector(double *v, int nl, int nh);
  96.     extern double *Ft_dvector(int nl, int nh);
  97.  
  98.     if (n != SplData) {
  99.         if (Xs)
  100.             Ft_free_dvector(Xs, 1, n);
  101.         if ((Xs = Ft_dvector(1, n)) == NULL)
  102.             return(ERRR);
  103.         if (Ys)
  104.             Ft_free_dvector(Ys, 1, n);
  105.         if ((Ys = Ft_dvector(1, n)) == NULL)
  106.             return(ERRR);
  107.         if (Y2s)
  108.             Ft_free_dvector(Y2s, 1, n);
  109.         if ((Y2s = Ft_dvector(1, n)) == NULL)
  110.             return(ERRR);
  111.     }
  112.     for (i=1;i<=n;i++) {
  113.         Xs[i] = x[i];
  114.         Ys[i] = y[i];
  115.     }
  116.     return(0);
  117. }
  118.  
  119.